In [1]:
import astropy.io.fits as fits
import numpy as np
from glob import glob
from matplotlib import pyplot as plt
import matplotlib.colors as colors
from astropy.visualization import simple_norm
import ccdproc as ccd
from astropy import visualization as viz
from astropy.visualization import (imshow_norm, MinMaxInterval, SqrtStretch, ZScaleInterval)
from astropy.coordinates import SkyCoord
from photutils.aperture import ApertureStats
from astropy.stats import SigmaClip
from photutils.aperture import aperture_photometry, CircularAperture, CircularAnnulus
from photutils.detection import DAOStarFinder
from astropy.stats import mad_std
from astropy.wcs import WCS
from astropy.table import Table
from astropy.io import ascii
from astropy import units as u
from datetime import datetime
from astropy.time import Time
In [2]:
# Sunset is 7:09 PM 
# By at observatory at 6:39 PM to take sky flats
# Ideally be done by 10 PM 
In [3]:
data = ascii.read('landoltstars.txt')  
In [4]:
star_pos = SkyCoord('16h41m41.24s', '36:27:35.5', unit=(u.hourangle, u.deg))
star_pos
Out[4]:
<SkyCoord (ICRS): (ra, dec) in deg
    (250.42183333, 36.45986111)>
In [5]:
data[0]
Out[5]:
Row index=0
StarDesignationRA(2000)Dec(2000)VB-VU-BV-RR-IV-I
str11str8str9float64float64float64float64float64float64
TPHE A00:30:09-46:31:2214.6510.7930.380.4350.4050.841
In [6]:
len(data)
Out[6]:
526
In [7]:
st = SkyCoord(data[0]['RA(2000)'], data[0]['Dec(2000)'], unit=(u.hourangle, u.deg))
st
Out[7]:
<SkyCoord (ICRS): (ra, dec) in deg
    (7.5375, -46.52277778)>
In [8]:
st.separation(star_pos).hour
Out[8]:
8.874273021712979
In [9]:
lands = []
for star in data:
    ra = star['RA(2000)']
    dec = star['Dec(2000)']
    coord = SkyCoord(ra, dec, unit=(u.hourangle, u.deg))
    #dist.append(coord.separation(star_pos))
    if star['V'] < 9.5 and coord.separation(star_pos).degree < 60:
        print(star)
        print('Separation:',coord.separation(star_pos).degree)
        lands.append(star)
StarDesignation RA(2000) Dec(2000)   V    B-V   U-B   V-R   R-I   V-I 
--------------- -------- --------- ----- ----- ----- ----- ----- -----
        109 231 17:45:20 -00:26:13 9.332 1.462 1.593 0.785 0.704 1.492
Separation: 39.74439851862938
StarDesignation RA(2000) Dec(2000)   V    B-V   U-B   V-R   R-I   V-I 
--------------- -------- --------- ----- ----- ----- ----- ----- -----
        111 773 19:37:17 +00:10:38 8.963 0.206 -0.21 0.119 0.144 0.262
Separation: 54.45338332240191
In [10]:
lands
Out[10]:
[<Row index=388>
 StarDesignation RA(2000) Dec(2000)    V      B-V     U-B     V-R     R-I     V-I  
      str11        str8      str9   float64 float64 float64 float64 float64 float64
 --------------- -------- --------- ------- ------- ------- ------- ------- -------
         109 231 17:45:20 -00:26:13   9.332   1.462   1.593   0.785   0.704   1.492,
 <Row index=431>
 StarDesignation RA(2000) Dec(2000)    V      B-V     U-B     V-R     R-I     V-I  
      str11        str8      str9   float64 float64 float64 float64 float64 float64
 --------------- -------- --------- ------- ------- ------- ------- ------- -------
         111 773 19:37:17 +00:10:38   8.963   0.206   -0.21   0.119   0.144   0.262]

Landolt Star Photometry¶

$V_{Landolt} = V_{Observed} + a (X-1) + b (X-1) (V_{Landolt}-I_{Landolt}) + ZP$

In [11]:
# For standard correction
# air mass
# airmass and color
# vega? 
In [ ]:
 
In [ ]:
 
In [ ]:
 

Data Reduction¶

In [12]:
directory = '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/*'
glob(directory)
Out[12]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/bias',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/dark',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/flat_I',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/flat_V',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/Observation Plan (10_6).docx',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/ObservationPlan.pdf',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/vega']
In [13]:
files_bias = glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/bias/*')
files_dark = glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/dark/*')
files_flatI =  glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/flat_I/*')
files_flatV =  glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/flat_V/*')
#files_m2 =  glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/*')
files_m2_I = glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/*')
files_m2_V = glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Vband/*')
files_vega_I =  glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/vega/*_I_*')
files_vega_V =  glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/vega/*_V_*')
In [14]:
def get_exptime(filename):
    return fits.getheader(filename)['EXPTIME']
def get_filter(filename):
    return fits.getheader(filename)['FILTER']
def get_obstime(filename):
    return Time(fits.getheader(filename)['DATE-OBS'])
def get_localtime(filename):
    return fits.getheader(filename)['NOTE'][11:]
In [15]:
names = ['Target', 'Exposure Time (s)', 'Filters', 'Number of Exposures', 'ObsTime (UTC)', 'Right Ascension', 'Declination']
filetype = ['Bias', 'Dark', 'I-Band Flat', 'V-Band Flat', 'Messier 2 I-Band', 'Messier 2 V-Band', 'Vega I-Band', 'Vega V-Band']
filelist = [files_bias, files_dark, files_flatI, files_flatV, files_m2_I, files_m2_V, files_vega_I, files_vega_V]
exptime = []
filters = []
numfiles = []
obstime = []
localtime = []
coord_m2 = SkyCoord('21h33m27s', '-00d49m23.7s', frame='fk5')
coord_vega = SkyCoord('18h36m56s', '+38d47m1s', frame='fk5')

coords_ra  = ['-', '-', '-', '-', '21h33m27s', '21h33m27s', '18h36m56s', '18h36m56s']
coords_dec = ['-', '-', '-', '-', '-00$^{\circ}$49\'23.7\"', '-00$^{\circ}$49\'23.7\"', '+38$^{\circ}$47\'1\"', '+38$^{\circ}$47\'1\"']

for lis in filelist:
    exptime.append(get_exptime(lis[0]))
    filters.append(get_filter(lis[0]))
    numfiles.append(len(lis))
    obstime.append(get_obstime(lis[0]))
    localtime.append(get_localtime(lis[0]))
In [16]:
Time(fits.getheader(files_m2_V[0])['DATE-OBS'])
Out[16]:
<Time object: scale='utc' format='isot' value=2022-10-07T01:51:05.000>
In [ ]:
 
In [17]:
fits.getheader(files_bias[0])['NOTE'][11:]
Out[17]:
'10/6/2022 at 18:57:46.000'
In [18]:
tblr = Table([filetype, exptime, filters, numfiles, obstime, coords_ra, coords_dec], names=names)
tblr.write('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/obstable.dat', format='latex', overwrite=True)
In [19]:
pwd
Out[19]:
'/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive'
In [20]:
fits.getheader(files_bias[0])
Out[20]:
SIMPLE  =                    T /                                                
BITPIX  =                   16 /                                                
NAXIS   =                    2 /                                                
NAXIS1  =                  765 /                                                
NAXIS2  =                  510 /                                                
OBJECT  = '        '           /                                                
TELESCOP= 'Celestron CPC 12'''                                                  
INSTRUME= 'SBIG ST-402'                                                         
OBSERVER= 'Daniel_Emma_Alex'                                                    
DATE-OBS= '2022-10-06T22:57:46.000' / GMT START OF EXPOSURE [WIN]               
BZERO   = +3.276800000000E+004 /                                                
BSCALE  = +1.000000000000E+000 /                                                
EXPTIME = +4.000000000000E-002 / EXPOSURE IN SECONDS                            
CCD-TEMP= -5.232156845990E+000 / CCD TEMP IN DEGREES C                          
XPIXSZ  = +9.000000000000E+000 / PIXEL WIDTH IN MICRONS                         
YPIXSZ  = +9.000000000000E+000 / PIXEL HEIGHT IN MICRONS                        
XBINNING=                    1 / HORIZONTAL BINNING FACTOR                      
YBINNING=                    1 / VERTICAL BINNING FACTOR                        
XORGSUBF=                    0 / SUB_FRAME ORIGIN X POS                         
YORGSUBF=                    0 / SUB-FRAME ORIGIN Y POS                         
EGAIN   = +1.550000000000E+000 / ELECTRONS PER ADU                              
FOCALLEN= +3.048000000000E+003 / FOCAL LENGTH IN MM                             
APTDIA  = +3.048000000000E+002 / APERTURE DIAMETER MM                           
APTAREA = +6.840554060800E+004 / APERTURE AREA IN SQ-MM                         
CBLACK  =                 1044 / BLACK ADU FOR DISPLAY                          
CWHITE  =                 1173 / WHITE ADU FOR DISPLAY                          
PEDESTAL=                 -100 / ADD TO ADU FOR 0-BASE                          
DATAMAX =                65535 / SATURATION LEVEL                               
SBSTDVER= 'SBFITSEXT Version 1.0' / SBIG FITS EXTENSIONS VER                    
SWACQUIR= 'WinOPS Ver 5.61 Build 0-NT' / DATA ACQ SOFTWARE                      
SWCREATE= 'SBIG Win CCDOPS Version 5.61 Build 0-NT'                             
FILTER  = 'Photometric I'      / OPTICAL FILTER NAME                            
SNAPSHOT=                    1 / NUMBER IMAGES COADDED                          
DATE    = '2022-10-06'         / GMT DATE WHEN THIS FILE CREATED                
RESMODE =                    0 / RESOLUTION MODE                                
EXPSTATE= '125     '           / EXPOSURE STATE (HEX)                           
RESPONSE= +3.000000000000E+003 / CCD RESPONSE FACTOR                            
NOTE    = 'Local time:10/6/2022 at 18:57:46.000'                                
In [21]:
ccd_shape = fits.getdata(files_bias[0]).shape

Bias¶

In [22]:
bias = np.zeros((len(files_bias), ccd_shape[0], ccd_shape[1]))
for bb in range(len(files_bias)):
    #print(files_bias[bb])
    bias[bb] = fits.getdata(files_bias[bb])
In [23]:
master_bias = np.median(bias, axis=0)
plt.imshow(master_bias)
Out[23]:
<matplotlib.image.AxesImage at 0x7f2e06cee490>
In [24]:
master_bias.mean()
Out[24]:
1088.0907497116493
In [ ]:
 

Darks¶

In [25]:
exptime_darkI = fits.getheader(files_m2_V[0])['EXPTIME']
exptime_darkI
Out[25]:
20.0
In [26]:
exptime_darkV = fits.getheader(files_dark[0])['EXPTIME']
exptime_darkV
Out[26]:
20.0
In [27]:
dark = np.zeros((len(files_dark), ccd_shape[0], ccd_shape[1]))
for bb in range(len(files_dark)):
    #print(files_bias[bb])
    dark[bb] = fits.getdata(files_dark[bb])
In [28]:
master_dark = np.median(dark-master_bias, axis=0)
plt.imshow(master_dark)
np.mean(master_dark)
Out[28]:
9.130502370882994

Flats¶

In [29]:
exptime_flatI = fits.getheader(files_flatI[0])['EXPTIME']
In [30]:
flatI = np.zeros((len(files_flatI), ccd_shape[0], ccd_shape[1]))
for bb in range(len(files_flatI)):
    #print(files_bias[bb])
    flatI[bb] = fits.getdata(files_flatI[bb]) - master_bias - exptime_flatI/exptime_darkI * master_dark
In [31]:
master_flatI = np.median(flatI, axis=0)
plt.imshow(master_flatI)
np.mean(master_flatI)
Out[31]:
30056.210125195445
In [32]:
master_flatI_norm = master_flatI/np.mean(master_flatI)
plt.imshow(master_flatI_norm)
np.mean(master_flatI_norm)
Out[32]:
0.9999999999999997
In [33]:
exptime_flatV = fits.getheader(files_flatV[0])['EXPTIME']
In [34]:
flatV = np.zeros((len(files_flatV), ccd_shape[0], ccd_shape[1]))
for bb in range(len(files_flatV)):
    #print(files_bias[bb])
    flatV[bb] = fits.getdata(files_flatV[bb]) - master_bias - exptime_flatV/exptime_darkV * master_dark
In [35]:
master_flatV = np.median(flatV, axis=0)-master_bias
plt.imshow(master_flatV)
np.mean(master_flatV)
Out[35]:
25671.547239965388
In [36]:
master_flatV_norm = master_flatV/np.mean(master_flatV)
plt.imshow(master_flatV_norm)
np.mean(master_flatV_norm)
Out[36]:
1.0000000000000002

A Note about the flats:

These flats were taken as sky flats before astronomical twilight. As part of the process of aiming the telescope at the North star to center it, the CCD was removed from the back of the telescope. Afterwards, it was placed back on, but imperfectly so that the dust on the optics was not in the same place. In the calibrated science images to follow, the "donuts" from dust particles on the optics are visible because of this.

Calibration¶

$[(flat-bias)-\frac{t_f}{t_d}(dark-bias)]_{ave}[\frac{(object-bias)-\frac{t_o}{t_d}(dark-bias)}{(flat-bias)-\frac{t_f}{t_d}(dark-bias)}]_{pix}$

$[master\_flat-\frac{exptime\_flat}{exptime\_dark}(master\_dark)]_{ave}[\frac{(object-master\_bias) -\frac{exptime\_obj}{exptime\_dark}(master\_dark)}{(master\_flat)-\frac{exptime\_flat}{exptime\_dark}(master\_dark)}]_{pix}$

M2¶

In [37]:
files_m2_I
Out[37]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_001.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_002.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_003.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_004.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_001.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_002.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_003.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_004.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_005.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_006.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_007.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_008.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_009.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_010.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_011.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_012.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_013.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_014.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_015.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_016.FIT']
In [38]:
data_m2_I = np.zeros((len(files_m2_I), ccd_shape[0], ccd_shape[1]))

for ff in range(len(files_m2_I)):
    data_m2_I[ff] = (fits.getdata(files_m2_I[ff])-master_dark) / master_flatI_norm
In [39]:
pwd
Out[39]:
'/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive'
In [40]:
norm_ex = np.median(data_m2_I, axis=0)
imshow_norm(norm_ex, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
print('Very blurry because the telescope does not track well. Need to have something that accounts for the movement of the stars in every science frame.')
plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/blurry.pdf')
Very blurry because the telescope does not track well. Need to have something that accounts for the movement of the stars in every science frame.
In [41]:
im1 = data_m2_I[0]
imshow_norm(im1, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[41]:
(<matplotlib.image.AxesImage at 0x7f2df564f370>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2df56e06d0>)
In [42]:
bkg_sigma = mad_std(im1)  
daofind = DAOStarFinder(fwhm=4., threshold=3. * bkg_sigma)  
sources = daofind(im1)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
#print(sources)
In [43]:
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=5.)  
In [44]:
fig = plt.figure(figsize=(8,9))
ax = imshow_norm(im1, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot(color='blue', lw=1.5, alpha=0.5)
#annulus_aperture.plot(color='blue', lw=1.5, alpha=0.1)
print('')

In [45]:
star = (677, 243)#(597, 34)
In [46]:
#dist = np.sqrt((positions[:,0]-star[0])**2 + (positions[:,1]-star[1])**2)
#index = np.where(np.min(dist) == dist)
#pos = positions[index][0]
#star_pos.append(pos)
#ids.append(index)
In [47]:
#pos
In [48]:
from astropy.visualization import ImageNormalize
Image Combination¶
In [49]:
#imshow_norm(im1, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
#shifted = []
data_m2_I = np.zeros((len(files_m2_I), ccd_shape[0], ccd_shape[1]))

for ff in range(len(files_m2_I)):
    im = (fits.getdata(files_m2_I[ff])-master_dark) / master_flatI_norm
    bkg_sigma = mad_std(im)  
    daofind = DAOStarFinder(fwhm=4., threshold=3. * bkg_sigma)  
    sources = daofind(im)
    positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
    apertures = CircularAperture(positions, r=5.) 
    dist = np.sqrt((positions[:,0]-star[0])**2 + (positions[:,1]-star[1])**2)
    index = np.where(np.min(dist) == dist)
    pos = positions[index][0]
    aper = CircularAperture(positions, r=5.) 
    #aper[index].plot(color='blue', lw=1.5, alpha=0.5)
    
    shi = np.roll(im, (int(star[1]-pos[1]), int(star[0]-pos[0])), axis=(0,1))
    
    #fig, axs = plt.subplots(1, 2, figsize=(9,9))
    #axs[0].imshow(im, norm=ImageNormalize(data=im, interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
    #axs[0].scatter(pos[0], pos[1], marker='.', color='blue', alpha=0.5)
    #axs[0].scatter(star[0], star[1], marker='.', color='red', alpha=0.5)
    #axs[1].imshow(shi, norm=ImageNormalize(data=shi, interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
    #axs[1].scatter(pos[0], pos[1], marker='.', color='blue', alpha=0.5)
    #axs[1].scatter(star[0], star[1], marker='.', color='red', alpha=0.5)
    
    data_m2_I[ff] = shi
    #print(im)
    
    #np.roll(im, pos-star)
In [50]:
#star_V = (594, 290)
star_V = (555, 70)
#star_V = (355, 177)
#star_V = (600, 21)
#star_V = (599, 479)
#star_V = (682, 230)
In [51]:
#imshow_norm(fits.getdata(files_m2_V[0]), origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
#shifted = []
data_m2_V = np.zeros((len(files_m2_V), ccd_shape[0], ccd_shape[1]))

for ff in range(len(files_m2_V)):
    im = (fits.getdata(files_m2_V[ff])-master_dark) / master_flatV_norm
    bkg_sigma = mad_std(im)  
    daofind = DAOStarFinder(fwhm=4., threshold=2.5 * bkg_sigma)  
    sources = daofind(im)
    positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
    apertures = CircularAperture(positions, r=5.) 
    dist = np.sqrt((positions[:,0]-star_V[0])**2 + (positions[:,1]-star_V[1])**2)
    index = np.where(np.min(dist) == dist)
    pos = positions[index][0]
    aper = CircularAperture(positions, r=5.) 
    #aper[index].plot(color='blue', lw=1.5, alpha=0.5)
    
    print(np.sqrt((pos[0]-star_V[0])**2 + (pos[1]-star_V[1])**2))
    shi = np.roll(im, (int(star_V[1]-pos[1]), int(star_V[0]-pos[0])), axis=(0,1))
    
    #fig, axs = plt.subplots(1, 2, figsize=(9,9))
    #axs[0].imshow(im, norm=ImageNormalize(data=im, interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
    #axs[0].scatter(pos[0], pos[1], marker='.', color='blue', alpha=0.5)
    #axs[0].scatter(star_V[0], star_V[1], marker='.', color='red', alpha=0.5)
    #axs[1].imshow(shi, norm=ImageNormalize(data=shi, interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
    #axs[1].scatter(pos[0], pos[1], marker='.', color='blue', alpha=0.5)
    #axs[1].scatter(star_V[0], star_V[1], marker='.', color='red', alpha=0.5)
    
    if np.sqrt((pos[0]-star_V[0])**2 + (pos[1]-star_V[1])**2) > 30.:
        print('Selected star too far away')
    else: 
        data_m2_V[ff] = shi
    #print(im)
    
    #np.roll(im, pos-star)
26.53738398585733
21.83304490504877
21.953036079431303
46.23193614817681
Selected star too far away
12.23206935111877
10.70735201945291
9.572564077172094
9.091323261693075
8.902041990812997
43.06073330763313
Selected star too far away
1.3256166289841604
5.831531623950198
11.455285389002357
5.18858963701437
5.318678335272187
5.538186868950662
14.841000883137074
14.090627098947493
11.151455903061647
9.690312857984582
In [52]:
image_combined_I = np.median(data_m2_I, axis=0)
fig = plt.figure(figsize=(8,9))
ax = imshow_norm(image_combined_I, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
In [53]:
fig = plt.figure(figsize=(8,9))
ax = imshow_norm(data_m2_I[0], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
In [54]:
image_combined_V = np.median(data_m2_V, axis=0)
fig = plt.figure(figsize=(8,9))
ax = imshow_norm(image_combined_V, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
In [55]:
fig = plt.figure(figsize=(8,9))
ax = imshow_norm(data_m2_V[0], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
In [56]:
bkg_sigma = mad_std(image_combined_I)  
daofind = DAOStarFinder(fwhm=4., threshold=3. * bkg_sigma)  
sources = daofind(image_combined_I)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=5.)  
annuli = CircularAnnulus(positions, 5., 8.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(image_combined_I, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot(color='blue', lw=1.5, alpha=0.2)
annuli.plot(color='blue', lw=1.5, alpha=0.2)

print('')

In [57]:
aperstats = ApertureStats(image_combined_I, annuli)
bkg_mean = aperstats.mean
phot_table_I = aperture_photometry(image_combined_I, apertures)
aperture_area = apertures.area_overlap(image_combined_I)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_I['aperture_sum'] - total_bkg
phot_table_I['total_bkg'] = total_bkg
phot_table_I['aperture_sum_bkgsub'] = phot_bkgsub
In [58]:
bkg_sigma = mad_std(image_combined_V)  
daofind = DAOStarFinder(fwhm=4., threshold=3. * bkg_sigma)  
sources = daofind(image_combined_V)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=5.)  
annuli = CircularAnnulus(positions, 5., 8.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(image_combined_V, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot(color='blue', lw=1.5, alpha=0.2)
annuli.plot(color='blue', lw=1.5, alpha=0.2)

print('')

In [59]:
aperstats = ApertureStats(image_combined_V, annuli)
bkg_mean = aperstats.mean
phot_table_V = aperture_photometry(image_combined_V, apertures)
aperture_area = apertures.area_overlap(image_combined_V)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_V['aperture_sum'] - total_bkg
phot_table_V['total_bkg'] = total_bkg
phot_table_V['aperture_sum_bkgsub'] = phot_bkgsub
In [60]:
#color = phot_table_V['aperture_sum_bkgsub'] - phot_table_I['aperture_sum_bkgsub']
#magnitude = phot_table_V['aperture_sum_bkgsub']
In [61]:
x_VIshift = phot_table_V[4]['xcenter'] - phot_table_I[4]['xcenter']
y_VIshift = phot_table_V[4]['ycenter'] - phot_table_I[4]['ycenter']
In [62]:
x_VIshift
Out[62]:
$5.1888917 \; \mathrm{pix}$
In [63]:
y_VIshift
Out[63]:
$-12.35472 \; \mathrm{pix}$
In [64]:
int(star_V[1]-pos[1])
Out[64]:
-8
In [65]:
r = np.roll(image_combined_V, (-int(x_VIshift.value), -int(y_VIshift.value)), axis=(1,0))
plt.imshow(r)
Out[65]:
<matplotlib.image.AxesImage at 0x7f2df42fa7f0>
Photometry¶
In [66]:
bkg_sigma = mad_std(r)  
daofind = DAOStarFinder(fwhm=4., threshold=3. * bkg_sigma)  
sources = daofind(r)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=5.)  
annuli = CircularAnnulus(positions, 5., 8.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(r, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot(color='blue', lw=1.5, alpha=0.2)
annuli.plot(color='blue', lw=1.5, alpha=0.2)

print('')

In [67]:
aperstats = ApertureStats(image_combined_V, annuli)
bkg_mean = aperstats.mean
phot_table_V = aperture_photometry(image_combined_V, apertures)
aperture_area = apertures.area_overlap(image_combined_V)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_V['aperture_sum'] - total_bkg
phot_table_V['total_bkg'] = total_bkg
phot_table_V['aperture_sum_bkgsub'] = phot_bkgsub
In [68]:
nearest_star = np.zeros(len(phot_table_V))
for st in range(len(phot_table_V)):
    xV = phot_table_V[st]['xcenter']
    yV = phot_table_V[st]['ycenter']
    nearest = 999
    
    for ar in range(len(phot_table_I)):
        xI = phot_table_I[ar]['xcenter']
        yI = phot_table_I[ar]['ycenter']
        dist = np.sqrt((xV-xI)**2+(yV-yI)**2).value
        
        if dist < nearest:
            nearest_star[st] = int(ar)
            nearest = dist
In [69]:
nearest_star
Out[69]:
array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,  15.,   9.,
        15.,  10.,  11.,  12.,  13.,  15.,  16.,  17.,  20.,  21.,  22.,
        25.,  24.,  25.,  26.,  27.,  28.,  30.,  29.,  31.,  32.,  30.,
        31.,  33.,  34.,  35.,  36.,  37.,  38.,  39.,  42.,  41.,  44.,
        43.,  45.,  46.,  49.,  48.,  50.,  51.,  52.,  54.,  55.,  59.,
        57.,  60.,  58.,  61.,  62.,  63.,  64.,  65.,  66.,  74.,  69.,
        67.,  68.,  71.,  70.,  72.,  74.,  73.,  75.,  76.,  82.,  78.,
        77.,  81.,  80.,  83.,  84.,  85.,  86.,  87., 105.,  89.,  90.,
        91.,  92.,  94.,  93.,  95.,  96.,  97.,  96.,  99.,  98., 103.,
       102., 104., 105., 109., 106., 111., 107., 108., 110., 112., 114.,
       113., 116., 115., 118., 120., 117., 121., 119., 122., 123., 124.,
       125., 115., 126., 126., 127., 127., 128., 129., 130., 131., 133.,
       132., 136., 134., 135., 138., 137., 139., 141., 140., 142., 151.,
       143., 144., 146., 147., 148., 144., 149., 150., 151., 152., 153.,
       154., 155., 157., 156., 158., 159., 160., 162., 161., 163., 164.,
       165., 167., 166., 168., 169., 170., 171., 172., 173., 175., 174.,
       177., 179., 180., 168., 182., 184., 183., 185., 186., 189., 188.,
       190., 191., 192., 194., 195., 196., 197., 198., 200., 201., 202.,
       203., 204., 205., 206., 207., 208.])
In [70]:
#color = phot_table_V['aperture_sum_bkgsub'] - phot_table_I['aperture_sum_bkgsub']
#magnitude = phot_table_V['aperture_sum_bkgsub']
In [71]:
nearest_star.size
Out[71]:
204
In [72]:
mag = []
color = []
for ss in range(nearest_star.size):
    xs = phot_table_V[ss]['xcenter'] - phot_table_I[int(nearest_star[ss])]['xcenter']
    ys = phot_table_V[ss]['ycenter'] - phot_table_I[int(nearest_star[ss])]['ycenter']
    dist = np.sqrt(xs**2 + ys**2).value
    #print(dist)
    if dist < 1.0:
        mag.append(phot_table_V[ss]['aperture_sum_bkgsub'])
        color.append(phot_table_V[ss]['aperture_sum_bkgsub'] - phot_table_I[int(nearest_star[ss])]['aperture_sum_bkgsub'])
In [73]:
plt.scatter(color, mag)
Out[73]:
<matplotlib.collections.PathCollection at 0x7f2df4123e50>
In [74]:
image_combined_I.shape
Out[74]:
(510, 765)
In [75]:
fig = plt.figure(figsize=(8,9))
ax = imshow_norm(image_combined_I, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot()
CircularAperture((295, 310), r=125.).plot(color='blue', alpha=0.5)
Out[75]:
(<matplotlib.patches.Circle at 0x7f2df406aa30>,)
In [76]:
mask = np.ones(image_combined_I.shape)
mask_bool = np.ones(image_combined_I.shape)
center = (295, 310)
radius = 125.
for x in range(mask.shape[1]):
    for y in range(mask.shape[0]):
        dist = np.sqrt((center[0] - x)**2+(center[1] - y)**2)
        if dist <= radius:
            mask[y,x] = 0#np.nan
            mask_bool[y,x] = np.nan
        dist2 = np.sqrt((mask.shape[1]/2 - x)**2+(mask.shape[0]/2 - y)**2)
        if dist2 >= 450:
            mask[y,x] = 0#np.nan
            mask_bool[y,x] = np.nan
        #print(dist)

plt.imshow(mask)
Out[76]:
<matplotlib.image.AxesImage at 0x7f2dde8fc100>
In [77]:
fig = plt.figure(figsize=(8,9))
ax = imshow_norm(r, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
CircularAperture((295, 310), r=125.).plot(color='blue', alpha=0.5)
CircularAperture((mask.shape[1]/2, mask.shape[0]/2), r=400.).plot(color='blue', alpha=0.5)
Out[77]:
(<matplotlib.patches.Circle at 0x7f2dde8e3310>,)
In [78]:
im_final_V = r*mask
im_final_I = image_combined_I*mask
im_final_I
Out[78]:
array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])
In [79]:
bkg_sigma = mad_std(im_final_V)  
daofind = DAOStarFinder(fwhm=4., threshold=3.0 * bkg_sigma)  
sources = daofind(im_final_V)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=5.)  
annuli = CircularAnnulus(positions, 5., 8.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(im_final_V*mask_bool, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot(color='blue', lw=1.5, alpha=0.2)
annuli.plot(color='blue', lw=1.5, alpha=0.2)

aperstats = ApertureStats(im_final_V, annuli)
bkg_mean = aperstats.mean
phot_table_V = aperture_photometry(im_final_V, apertures)
aperture_area = apertures.area_overlap(im_final_V)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_V['aperture_sum'] - total_bkg
phot_table_V['total_bkg'] = total_bkg
phot_table_V['aperture_sum_bkgsub'] = phot_bkgsub

print('')
plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/photometry.pdf', bbox_inches='tight')

In [80]:
bkg_sigma = mad_std(im_final_I)  
daofind = DAOStarFinder(fwhm=4., threshold=3.0 * bkg_sigma)  
sources = daofind(im_final_I)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=5.)  
annuli = CircularAnnulus(positions, 5., 8.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(im_final_I*mask_bool, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot(color='blue', lw=1.5, alpha=0.2)
annuli.plot(color='blue', lw=1.5, alpha=0.2)

aperstats = ApertureStats(im_final_I, annuli)
bkg_mean = aperstats.mean
phot_table_I = aperture_photometry(im_final_I, apertures)
aperture_area = apertures.area_overlap(im_final_I)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_I['aperture_sum'] - total_bkg
phot_table_I['total_bkg'] = total_bkg
phot_table_I['aperture_sum_bkgsub'] = phot_bkgsub

print('')

In [81]:
nearest_star = np.zeros(len(phot_table_V))
for st in range(len(phot_table_V)):
    xV = phot_table_V[st]['xcenter']
    yV = phot_table_V[st]['ycenter']
    nearest = 999
    
    for ar in range(len(phot_table_I)):
        xI = phot_table_I[ar]['xcenter']
        yI = phot_table_I[ar]['ycenter']
        dist = np.sqrt((xV-xI)**2+(yV-yI)**2).value
        
        if dist < nearest:
            nearest_star[st] = int(ar)
            nearest = dist
In [82]:
mag = []
color = []
for ss in range(nearest_star.size):
    xs = phot_table_V[ss]['xcenter'] - phot_table_I[int(nearest_star[ss])]['xcenter']
    ys = phot_table_V[ss]['ycenter'] - phot_table_I[int(nearest_star[ss])]['ycenter']
    dist = np.sqrt(xs**2 + ys**2).value
    #print(dist)
    if dist < 1.0:
        mag.append(phot_table_V[ss]['aperture_sum_bkgsub'])
        color.append(phot_table_V[ss]['aperture_sum_bkgsub'] - phot_table_I[int(nearest_star[ss])]['aperture_sum_bkgsub'])
In [83]:
plt.scatter(color, mag)
plt.gca().invert_yaxis()
In [84]:
fig = plt.figure(figsize=(8,9))
ax = imshow_norm(r, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
# Final V-band
plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/messier2V.pdf', bbox_inches='tight')
In [85]:
fig = plt.figure(figsize=(8,9))
ax = imshow_norm(image_combined_I, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
# Fianl I-Band
plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/messier2I.pdf', bbox_inches='tight')

Vega¶

$[master\_flat-\frac{exptime\_flat}{exptime\_dark}(master\_dark)]_{ave}[\frac{(object-master\_bias) -\frac{exptime\_obj}{exptime\_dark}(master\_dark)}{(master\_flat)-\frac{exptime\_flat}{exptime\_dark}(master\_dark)}]_{pix}$

Image Combination¶
In [86]:
exptime_vegaI = fits.getheader(files_vega_I[0])['EXPTIME']
In [87]:
exptime_darkI
Out[87]:
20.0
In [88]:
data_vega_I = np.zeros((len(files_vega_I), ccd_shape[0], ccd_shape[1]))

for ff in range(len(files_vega_I)):
    data_vega_I[ff] = (fits.getdata(files_vega_I[ff])-(exptime_vegaI/exptime_darkI)*master_dark) / master_flatI_norm
In [89]:
vega_final_I = np.median(data_vega_I, axis=0)
imshow_norm(vega_final_I, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[89]:
(<matplotlib.image.AxesImage at 0x7f2df43e4940>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ddeac5eb0>)
In [90]:
exptime_vegaV = fits.getheader(files_vega_V[0])['EXPTIME']
In [91]:
data_vega_V = np.zeros((len(files_vega_V), ccd_shape[0], ccd_shape[1]))

for ff in range(len(files_vega_V)):
    data_vega_V[ff] = (fits.getdata(files_vega_V[ff])-(exptime_vegaV/exptime_darkV)*master_dark) / master_flatV_norm
In [92]:
vega_final_V = np.median(data_vega_V, axis=0)
imshow_norm(vega_final_V, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[92]:
(<matplotlib.image.AxesImage at 0x7f2df42a4d00>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2df4396fd0>)

Note: Vega is definitely over exposed

Photometry¶
In [93]:
bkg_sigma = mad_std(vega_final_I)  
daofind = DAOStarFinder(fwhm=15., threshold=100.0 * bkg_sigma)  
sources = daofind(vega_final_I)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=60.)  
annuli = CircularAnnulus(positions, 5., 8.)

ap_I = CircularAperture((498, 197), r=60.)
an_I = CircularAnnulus((498, 197), 70., 120.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(vega_final_I, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')

ap_I.plot(color='red', lw=1.5, alpha=0.2)
an_I.plot(color='red', lw=1.5, alpha=0.2)

#apertures.plot(color='blue', lw=1.5, alpha=0.2)
#annuli.plot(color='blue', lw=1.5, alpha=0.2)

print('')
plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/vegaI.pdf', bbox_inches='tight')

In [94]:
aperstats = ApertureStats(vega_final_I, an_I)
bkg_mean = aperstats.mean
phot_table_vegaI = aperture_photometry(vega_final_I, ap_I)
aperture_area = ap_I.area_overlap(vega_final_I)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_vegaI['aperture_sum'] - total_bkg
phot_table_vegaI['total_bkg'] = total_bkg
phot_table_vegaI['aperture_sum_bkgsub'] = phot_bkgsub

phot_table_vegaI
Out[94]:
QTable length=1
idxcenterycenteraperture_sumtotal_bkgaperture_sum_bkgsub
pixpix
int64float64float64float64float64float64
1498.0197.059047290.5342532412600461.47632765846446829.05792558
In [ ]:
 
In [95]:
bkg_sigma = mad_std(vega_final_V)  
daofind = DAOStarFinder(fwhm=16., threshold=100.0 * bkg_sigma)  
sources = daofind(vega_final_V)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=80.)  
annuli = CircularAnnulus(positions, 5., 8.)

ap_V = CircularAperture((469, 185), r=75.)
an_V = CircularAnnulus((469, 185), 100., 150.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(vega_final_V, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')

ap_V.plot(color='red', lw=1.5, alpha=0.2)
an_V.plot(color='red', lw=1.5, alpha=0.2)

#apertures.plot(color='blue', lw=1.5, alpha=0.2)
#annuli.plot(color='blue', lw=1.5, alpha=0.2)

print('')
plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/vegaV.pdf', bbox_inches='tight')

In [96]:
aperstats = ApertureStats(vega_final_V, an_V)
bkg_mean = aperstats.mean
phot_table_vegaV = aperture_photometry(vega_final_V, ap_V)
aperture_area = ap_V.area_overlap(vega_final_V)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_vegaV['aperture_sum'] - total_bkg
phot_table_vegaV['total_bkg'] = total_bkg
phot_table_vegaV['aperture_sum_bkgsub'] = phot_bkgsub

phot_table_vegaV
Out[96]:
QTable length=1
idxcenterycenteraperture_sumtotal_bkgaperture_sum_bkgsub
pixpix
int64float64float64float64float64float64
1469.0185.0115364750.2482222219942595.75966265895422154.48855956
In [97]:
phot_table_I
Out[97]:
QTable length=82
idxcenterycenteraperture_sumtotal_bkgaperture_sum_bkgsub
pixpix
int64float64float64float64float64float64
113.5540624491241811.84660995481611561986.2596530730653826.4551919042358159.804461168824
2751.42504616221451.851170779937263262041.9759183886353734.784572164738307.191346223895
3595.105232628340433.10370087945881100503.6093102833192804.31963971457699.289670568818
420.4199239654261158.7927925184922394673.8273802052692551.692449671782122.1349305334734
5256.862278139844365.62829102346362100167.175921722693532.44994334186634.725978380811
6337.987900399974477.3539535715546898702.9891333534393535.271153988435167.717979364999
7549.81648725534482.2000357425378298768.8839769422892853.51134682925915.372630113081
8230.9234083973787285.27939424051158101669.1238743027294218.681882451137450.441991851592
9333.474586226284890.52576834219148101319.2145367883393921.228386626727397.986150161611
10188.58014055544743121.733561987702496492.310133189193359.354259875563132.9558733135345
11512.7921102699793123.5917214558294797293.9326109540792840.212824662224453.719786291855
12349.57846761791666130.8166465939787113854.2375308289596632.5482751621217221.68925566683
13129.77302367955207131.243645256748199218.1019996338795290.84426289763927.2577367362683
14122.77009753307367135.189988216222497495.912116424895454.490094742322041.4220216824906
15278.1797807816903136.9606640195692496182.0940462297393668.099786375722513.9942598540074
..................
68204.03202978693398428.9174337003223699248.4169759978597790.328597025521458.0883789723302
69400.54870368739836428.943209378363798833.6354217654795432.283201607543401.3522201579326
70194.48150990383192433.9323301665131124338.46285111597100554.7607036125823783.702147503383
71195.09763830666935442.9130690818934799777.950489144100878.42210998737-1100.4716208433674
72137.49988682444481451.0578691384308101620.3573150909796701.389555810794918.9677592801745
73196.35965946954394457.9901744254795104313.7860582777696346.212431275857967.5736270019115
74280.9762331964428459.04722896464483125933.56636823025100638.4465978058925295.119770424368
75343.9128415509698460.0612440892582599254.243788084695425.64136571783828.602422366792
76133.54062267889867461.3818727820196102670.4946463463396426.590974357696243.903671988635
77622.8511800857611474.369425225922598102.7264624399294515.393909973143587.332552466789
78445.0018572895247483.227263620490599748.4186423931894133.628355144175614.7902872490085
79327.9762114437748488.519061703318102915.016255579595292.449203705667622.567051873833
80290.11404339062364489.8215848313834102748.7518605241895774.42100429216974.330856232074
81592.551732093303491.3400594929506107048.315400201795480.1053367809411568.21006342076
82186.57758077001452499.57231970824216112927.4800970787396145.6584477462716781.821649332458

Magnitude Calculations¶

m$_1 - $ m$_{ref} = -2.5 log_{10}(\frac{I_1}{I_{ref}})$

In [98]:
exptime_m2I = exptime_darkI
exptime_m2V = exptime_darkV
gain = fits.getheader(files_m2_I[0])['EGAIN'] # electrons per adu
In [99]:
intensity_m2V = phot_table_V['aperture_sum_bkgsub']*gain/exptime_m2V
intensity_m2I = phot_table_I['aperture_sum_bkgsub']*gain/exptime_m2I

intensity_vegaV = phot_table_vegaV['aperture_sum_bkgsub']*gain/exptime_vegaV
intensity_vegaI = phot_table_vegaI['aperture_sum_bkgsub']*gain/exptime_vegaI
In [100]:
mag_m2V = -2.5*np.log10(intensity_m2V/intensity_vegaV)
mag_m2I = -2.5*np.log10(intensity_m2I/intensity_vegaI)
/tmp/ipykernel_5060/3524271707.py:2: RuntimeWarning: invalid value encountered in log10
  mag_m2I = -2.5*np.log10(intensity_m2I/intensity_vegaI)
In [101]:
phot_table_V['magnitude'] = mag_m2V
phot_table_I['magnitude'] = mag_m2I
In [102]:
nearest_star = np.zeros(len(phot_table_V))
for st in range(len(phot_table_V)):
    xV = phot_table_V[st]['xcenter']
    yV = phot_table_V[st]['ycenter']
    nearest = 999
    
    for ar in range(len(phot_table_I)):
        xI = phot_table_I[ar]['xcenter']
        yI = phot_table_I[ar]['ycenter']
        dist = np.sqrt((xV-xI)**2+(yV-yI)**2).value
        
        if dist < nearest:
            nearest_star[st] = int(ar)
            nearest = dist
In [103]:
mag_V = []
mag_I = []
color = []
for ss in range(nearest_star.size):
    xs = phot_table_V[ss]['xcenter'] - phot_table_I[int(nearest_star[ss])]['xcenter']
    ys = phot_table_V[ss]['ycenter'] - phot_table_I[int(nearest_star[ss])]['ycenter']
    dist = np.sqrt(xs**2 + ys**2).value
    #print(dist)
    if dist < 1.0:
        mag_V.append(phot_table_V[ss]['magnitude'])
        mag_I.append(phot_table_I[int(nearest_star[ss])]['magnitude'])
        color.append(phot_table_V[ss]['magnitude'] - phot_table_I[int(nearest_star[ss])]['magnitude'])
In [104]:
fig, axs = plt.subplots(1, 2, figsize=(9,5))

axs[0].scatter(color, mag_I)
axs[0].invert_yaxis()
axs[0].set_xlabel('Color (V-I)')
axs[0].set_ylabel('I-band Magnitude')

axs[1].scatter(color, mag_V)
axs[1].invert_yaxis()
axs[1].set_xlabel('Color (V-I)')
axs[1].set_ylabel('V-band Magnitude')

plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/colormag.pdf')
In [ ]:
 
In [ ]: